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Abstract - We study numerically local models for the mechanical contact between two solids 
with rough surfaces. When the solids softly touch either through adhesion or by a small normal 
load L, contact only forms at isolated patches and fluids can pass through the interface. When the 
load surpasses a threshold value, L c , adjacent patches coalesce at a critical constriction, i.e., near 
points where the interfacial separation between the undeformed surfaces forms a saddle point. 
This process is continuous without adhesion and the interfacial separation near percolation is 
fully defined by scaling factors and the sign of L c — L. The scaling factors lead to a Reynolds 
flow resistance which diverges as (L c — L) -/3 with — 3.45. Contact merging and destruction 
near saddle points becomes discontinuous when either short-range adhesion or specific short-range 
repulsion are added to the hard-wall repulsion. These results imply that coalescence and break-up 
of contact patches can contribute to Coulomb friction and contact aging. 


Introduction. — When two nominally rough solids 
are pressed against each other, their surfaces tend to touch 
microscopically only at isolated points [l][2]. Simple mod¬ 
els assume that real contact between two rough surfaces 
can be decomposed into single-asperity, Hertzian-like con¬ 
tacts. However, both experiment [3|4] and large-scale sim¬ 
ulations [5j|9 have revealed that the majority fraction of 
contiguous contact patches appear to be fractal and to re¬ 
sult from many (formerly) single-asperity contacts. While 
the latter are well studied, less is known about the way 
in which they merge [To] 13 . in particular when adhesion 
is present. For example, it is unclear if adhesive contact- 
patch coalescence and break-up happen discontinuously. 
Contact-patch coalescence also plays a role for seals: it 
has been argued that their leakage rate is determined to 
a significant degree from the topology of the last criti¬ 
cal constriction II], i.e., the neighborhood of the point, 
which, upon increasing load or adhesion, is the first point 
to interrupt a percolating non-contact path. 

In this work, we study local models for the merging of 
contact patches, or, in the context of seals, the contact 
mechanics of critical constrictions. The local gap topog¬ 
raphy of a critical constriction is such that the interfacial 
separation of the undeformed surfaces is (close to) a sad¬ 
dle point. More precisely, near the percolation point, the 
gap is very small. Parallel to a just-blocked fluid channel, 


the gap opens, while in the orthogonal, in-plane direction, 
the gap is closed, because the interfacial separation of the 
undeformed surfaces decreases in that direction. Thus, 
studying the contact mechanics of saddle points entails 
the analysis of critical constrictions and that of contact 
patch merging. Here, we investigate not only the con¬ 
tact mechanics of (near-) critical constrictions, including 
a scaling analysis of the gap topography near the percola¬ 
tion point, but also address the pertaining Reynolds flow 
with and without adhesion. 

To simulate the contact mechanics of an isolated sad¬ 
dle point, we reinvestigate rather simple models fio][TT][T3l 
for the (combined) surface roughness and for the interac¬ 
tions between the surfaces. Yet, in addition to the hard- 
wall repulsion commonly assumed infinitely short-ranged 
in continuum mechanics, we also study the effects of 
finite-range attraction [l5] and repulsion 16 between the 
surfaces. Our motivation for additional repulsion stems 
from recent simulations, in which appropriately designed 
cohesive-zone models qualitatively reproduced the surface 
interactions mediated by a strongly wetting fluid inducing 
an effective negative surface energy of finite range fl6] . 
The length-scale of these additional surface interactions is 
nevertheless short enough for interfacial interactions not 
to act far away from a contact line, i.e., our Tabor param¬ 
eters 15] are >> 1. 
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Models and Methods. — We assume linear elastic¬ 
ity and the small-slope approximation so that roughness 
can be mapped to a rigid substrate and the elastic compli¬ 
ance to a flat counter body. The effective contact modulus 
is used to define the unit of pressure, i.e., E* = 1. Elas¬ 
ticity is treated with Green’s function molecular dynam¬ 
ics (GFMD) [l7 and the continuum version of the stress- 
displacement relation in Fourier space, cf(q) = qE*u(q)/2, 
where q is a wave vector and q its magnitude. Simulations 
are run in a force-controlled fashion. After the external 
load has changed by a small amount, all degrees of free¬ 
dom are relaxed until convergence is attained. 

Three models for the (combined) height profiles are em¬ 
ployed. Each profile has roughness only at a single wave¬ 
length A, a root-mean-square height fluctuation of ho, 
and a root-mean-square gradient of 27rho/A, which should 
be small enough for the small-slope approximation to be 
valid. In other words, all three lattices yield identical 
angle-averaged spectra at non-zero wave numbers. First, 
we consider a square lattice 
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for which the simulation cell dimensions along x and y 
direction are chosen to coincide with the wavelength A of 
the height undulation, that is L x = L y = A. The two 
remaining lattices are triangular and hexagonal: 
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for which the dimensions of the periodically repeated sim¬ 
ulation cells are L x = 4A/\/3 and L y = 2A. This shape is 
relatively close to a square so that discretization correc¬ 


tions 18 are almost isotropic in our GFMD calculations, 


which decomposes the surfaces into 2 n x 2 n elements. 

In each model, the height offset is set to yield a min¬ 
imum height of zero. Other key geometrical attributes 
are: Maximum height (y^27/2 ho, 4 ho, ^27/2 ho), mean 
trough depth (\/6 ho, 2 ho, y/2/3 ho), skewness (—^2/3, 
0, 2/3), and excess kurtosis (-1/2, -3/4, -1/2), each time 

in the order hexagonal, square, and triangular lattice. 

Two surfaces interact with a hard-wall constraint, i.e., 
they are not allowed to overlap. In addition we assume 
cohesive zone models for the finite-range interactions be¬ 
tween the surfaces. It is 7 = —70 exp{— g/ p} for attractive 
forces, where g = g(x,y) is the local gap or interfacial sep¬ 
aration. Typical values for both 70 and p are small num¬ 
bers in the appropriate unit system, e.g., 70 = 0.2 F*hg/A 
and p = 0.04 ho- To give one possible realization: If 
E* were 1 GPa, ho = 10 nm, and A = 1 p m, then 
70 = 20 mN/m, which reflects the order of magnitude for 



Fig. 1: Elastic displacements at zero load for the hexago¬ 
nal, square, and triangular substrate lattice, from left to right. 
Dark colors represent small displacements, that is, points of 
contact. These are enveloped by yellow halos indicative of lo¬ 
calized, adhesive necks. Thin black lines connect substrate 
height maxima, while thin white lines connect substrate height 
minima or troughs, which are marked by white points. Grey 
stars reflect saddle points. Only selected lines between nearest- 
neighbor saddle points are drawn (thick dotted lines). 


the surface energy of flat, non-polar, and chemically pas¬ 
sivated surfaces. Finite-range repulsion is modeled with 
a different functional form than adhesion, namely with 
l(x,y) = —70 exp{— g(x, y) 2 /2p 2 }. The detailed rationale 
for these functional forms is given elsewhere 16 . In brief: 


the exponential adhesive model mimics the large Tabor 
number limit quite efficiently while allowing one to deter¬ 
mine contact area in an unambiguous manner. The Gaus¬ 
sian repulsive interactions can simulate the effect of a last 
layer, which needs to be squeezed out before two surfaces 
touch. 

Figure [l] conveys an impression of our models. It shows 
the elastic displacement for the investigated symmetries 
in case of adhesion and zero external load. The dark ar¬ 
eas show the locations of small displacement and thereby 
the structure of initial contact patches. The halos around 
contact points reflect bulges near the contact line that are 
well known from single-asperity contacts with short-range 
adhesion. One can also recognize that the displacements 
are rather constant near the height minima. The latter he 
on lattices dual to the height maxima for the considered 
cases. For example, the height minima (or gap maxima) 
form a triangular lattice when the substrate heights lie on 
a hexagonal lattice, and vice versa. Saddle points form 
a kagome lattice in these two cases. For the square lat¬ 
tice, not only height maxima but also minima and saddle 
points each form a square lattice. 

Flow through the interface is described by Reynolds 
thin-film equation, which assumes a local conductivity 
that scales as the inverse third power of the gap. We use 
the hypre package 19 to solve the sparse linear system 


that the discretized Reynolds equation can be expressed 
as. We employ the solvers supplied with hypre using the 
CG (conjugate gradient), or GMRES (generalized mini¬ 
mal residual) methods 20 , each preconditioned using the 
PFMG method, which is a parallel semicoarsening multi¬ 
grid solver [21 . Our in-house code is MPI-parallelized and 
uses HDF5 for I/O. 
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Fig. 2: Relative contact area as a function of load for the 
square-lattice substrate and different finite-range interactions. 
Red (blue) lines refer to increasing (decreasing) external load 
as marked by arrows. Exponential adhesion (solid lines) and 
Gaussian repulsion (dashed lines) show the established hystere- 
ses at small contact area related to contact formation. Unless 
solids interact solely by hard-wall repulsion, additional hystere- 
ses appear at relative contact areas near 0.4. They reflect the 
contact-patch-coalescence and break-up instabilities. 


Results. — A central question addressed in this study 
is whether the merging of two contact patches happens 
continuously or discontinuously. Figure [ 2 ] confirms previ¬ 
ous findings 10,13 that percolation is continuous when 
surfaces interact solely with hard-wall repulsion. How¬ 
ever, once finite-range repulsion (with the functional form 
defined in the method section) or short-range attraction 
are added, contact patches and likewise fluid channels co¬ 
alesce and break up through instabilities. The behavior is 
qualitatively similar for the two remaining models, despite 
large quantitative differences. 


Contacts also coalesce discontinuously when driven en¬ 
tirely by adhesion. This is revealed in Figure [3] Despite 
similar conditions — identical adhesion and surface spec¬ 
tra — the quantitative differences between the three mod¬ 
els are clearly borne out: The hexagonal (triangular) lat¬ 
tice has by far the largest (smallest) propensity to form a 
coalesced or percolated contact area. This is because its 
peaks are rather blunt or obtuse (pointed or acute) and 
the ridges connecting the peaks have small (large) cur¬ 
vature. However, the hexagonal (triangular) lattice has 
the smallest tendency to go into full contact, because its 
troughs are rather deep (shallow) and the curvature nor¬ 
mal to the ridges are large (small) in magnitude. In the 
triangular lattice, the propensity to form full contact is 
even so large that partial, coalesced contact is not even 
metastable at finite adhesion and zero load. The square 
model is somewhere in between the two extreme cases. 


While adhesion-free/load-driven and adhesion- 
driven/load-free percolation are qualitatively different in 
that the former is continuous and the latter is discon¬ 
tinuous, they happen at similar relative contact areas. 


Fig. 3: Relative contact area as a function of surface energy 
70 for the three investigated models. The instabilities at lower 
values of 70 are related to percolation transitions. The single¬ 
wavelength, triangular substrate circumvents the percolated, 
partial contact regime by transitioning directly from isolated 
contact patches to full contact. 


The precise values in the purely load-driven case with 
only hard-wall repulsion are a p = 0.17826(11) (triangular 
lattice), 0.40185(6) (square lattice), and 0.67323(1) 
(hexagonal lattice). Interestingly, the average of these 
numbers (0.418) is close to the percolation threshold 
identified for self-affine, isotropic, randomly-rough and 
non-adhesive bodies (0.425) [22]. In the adhesion-driven 
case, the percolation of contact is triggered at relative 
contacts slightly below the given values of a p and ends 
at values slightly above a p . The reverse transitions from 
partial but coalesced to isolated contact patches also 
occur near the respective values for a p , though shifted to 
slightly smaller values. 

While only peripheral to this work, we briefly discuss the 
transition to full contact. In the purely force-driven case, 
full contact occurs above critical mean pressures of pf c = 
7rE*t/\ , where t is the mean trough depth (stated in the 
model section). This relation is readily obtained from the 
single-wavelength, full-contact solution 110]. (In the orig¬ 
inal work, an additional factor of y/2 was included, as the 
wavelength was defined to be \/y/2.) The purely adhesion- 
driven percolation roughly follows the force-driven perco¬ 
lation, e.g., the ratio 7 f c / 7 P is largest for the hexagonal 
(135) and smallest for the triangular ( 1 ) model. Pertinent 
ratios for the critical forces in the absence of adhesion are 
31 (hexagonal) and 1.4 (triangular). Here and in the fol¬ 
lowing, we do not address the reverse transition from full 
to partial contact for adhesive surfaces, because it would 
correspond to a Griffith-like fracture problem lacking an 
initial crack. As such, for our and related adhesive laws, 
full contact becomes unstable at tensile pressure < 70 /p. 
The other adhesion-related results presented here barely 
depend on the precise choice of p as long as p <C Hq. 

The above results reveal a quite significant influence of 
local structure on prefactors. Some of the differences can 
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Fig. 4: Gaps (in color, where bright colors indicate large gaps 
and dark color small gaps) and contact (white) near the perco¬ 
lation threshold for the three investigated geometries. Systems 
are arranged in the same way as in Figure [l] 

be rationalized from the contact and gap geometries near 
the percolation threshold, which are depicted in Figure [4] 
for simple hard-wall repulsion. In coarse-grained repre¬ 
sentations the impression is conveyed that, depending on 
geometry, opening angles can range from acute to obtuse. 

Locally, however, all investigated saddle points are sim¬ 
ilar in that they sit at points of inversion symmetry, the 
curvatures of the substrate heights parallel to the fluid- 
channel (/ijj) are negative and those in orthogonal, in¬ 
plane direction (h f [) are positive. One can therefore ex¬ 
pect that the solutions for gap and displacement are (lo¬ 
cally) similar for the different symmetries, except for pre- 
and scaling factors. In addition, percolation is continuous 
in the adhesion-free case, thereby constituting a critical 
point. Thus, there should exist (quasi-) universal func¬ 
tions for the various fields in the vicinity of the critical 
point, e.g., for the gap, on which we focus exemplarily. 
Following the ideas of the scaling hypothesis, the gap 
should then satisfy 

g(l,x,y) = \l\<g ±( 4 ) 

where l = (L — L c )/L c is a reduced load, g±(. . .) denotes 
two master functions depending on the sign of /, while x> 
v, and C are universal scaling exponents. To superimpose 
the gap functions, not only for different values of l but 
also for different models — or more generally speaking for 
different curvature ratios 77 = —h'^/h'[_ characterizing the 
geometry of the saddle point — appropriate unit systems 
for x, 7 /, and 2 need to be defined. For example, for the 
square lattice, or 77 = 1 , a reasonable choice is to define the 
minimum width and minimum height of the open channel 
as 1 for l = —0.01 just like the length of the blocked chan¬ 
nel for l = 0.01. For other geometries, or 77 7 ^ 1, the proper 
definition of units then depends on 77 and also on the sign 
of /. 

Numerically, we see equation © to hold and demon¬ 
strate this in Figure |5j Specifically, in panel (a), contact 
lines obtained for different reduced loads and different ge¬ 
ometries (exemplarily triangular and square) collapse onto 
their master curve that solely depends on the sign of 1. 
Thus, all (near-) critical constrictions from Figure [4j lo- 



distance from saddle point 

Fig. 5: (a) Contact line shape for normal loads below (blue), 

at (green), and, above (red) the critical load, (b) Gap on 
the symmetry axis as a function of the distance from the sad¬ 
dle point. Different symbols represent different saddle-point 
geometries and different colors different reduced loads l , as in¬ 
dicated in the caption. Solid areas and lines are obtained from 
analytical approximations to the contact shape. At a given 
value of l , distances are rescaled according to equation ©• 
Units are chosen such that gap height, length, and width are 1 
for |/| = 0.01 in the square model. 

cally adopt the contact line shapes from Figure [5] when 
viewed with sufficient magnification. Panel (b) shows that 
the gap can also be collapsed onto universal functions. 

The numerical estimates for the exponents were ob¬ 
tained as follows: We first identified L c as accurately as 
possible. We then ran simulations for l = —0.002, —0.005, 
—0.01, and —0.02. By analyzing how the gap at the sad¬ 
dle point scales with /, the exponent ( was determined to 
be C = 1.2 ± 0.01. The widths of the open channels were 
determined from the same set of simulations. The con¬ 
tact line positions, r c , resulted from fitting the relation 
F(r || = 0,rj_) oc y/r± — r c to the forces exerted by the 
substrate onto the elastic body near the contact line. Here, 
the saddle point sat at r*| | =0 and r± = 0 in the local coor¬ 
dinate system. This fitting procedure allowed us to deter¬ 
mine the width of the open channel with sub-discretization 
resolution. All data was consistent with v = 0.45 db 0.01. 
To determine the length of blocked channels, simulations 
were run for positive reduced loads and the same absolute 
values of l as before. The pressure profile on the symme¬ 
try axis (7*||,rj_ = 0) followed a Hertzian pressure profile 
very accurately, which allowed us to determine the length 
of the contact line on the symmetry axis to high precision 
and ultimately to ascertain that x = 0.6±0.01. Since the 
pressure profile on the contact ridges appears to be iden¬ 
tical to that of a Hertzian contact, we believe that there 
exists an analytical solution, in which case, we would ex¬ 
pect the exponents to be simple rational numbers — as in 
the Hertzian contact problem. We therefore speculate that 
X = 3/5, v = 9/20, and ( = 6/5 are the exact exponents. 


P-4 
















Contact mechanics of saddle points 


L 



Fig. 6 : Top: Fluid flow through the interface as a function of 
load. Solid lines refer to 70 = 0 and dashed lines to 70 = 0.05. 
The latter have end points, which are indicated by symbols. 
Bottom left: Normalized fluid flow as a function of L c — L, 
in the non-adhesive case. Bottom right: Visualization of the 
Reynolds thin-film current density in the critical region. 


Knowledge of the three exponents and v allows 

one to deduce the power law with which Reynolds flow 
is suppressed as the load is increased toward L c . Since 
the fluid-pressure gradient is non-negligible only near the 
critical constriction, the resistance to fluid flow is propor¬ 
tional to the length of the channel. Moreover, it is in¬ 
versely proportional to the channel width and, according 
to the Reynolds equation, the third power of the channel 
height. Thus, the resistance scales as 

R(L ) oc (L c - L)~P (5) 


with 


P = K + v-x, 


(6) 


and a numerical value of 13 = 3.45. Figure [6] confirms these 
expectations for the investigated adhesion-free surfaces: 
Near the percolation threshold, all three system disappear 
with the same power law, even if prefactors may differ 
substantially. The exponent of /3 = 3.45 is reproduced 
quite accurately by the flow simulations. In the critical 
regime, the (scaled) current density always looks like the 
flow shown in the bottom right of Figure [6] 

The absolute flow through the contacts changes dra¬ 
matically with the saddle point geometry even for fixed 
height spectra, in parts, of course, because the percola¬ 
tion point is so sensitive to the geometry. In other words, 
the leakage rate is strongly dependent on the curvature 
ratio Tj or the skewness s of the height distribution, po¬ 
tentially even more so than previously recognized 23 , see 
also 24 . Specifically, the single-wavelength models have 
f) = (4,1,1/4) and s = (—y/2/3,0, y/2/3) for hexagonal, 
square, and triangular lattices, respectively. In addition, 
just as the contact area is discontinuous at the percolation 


transition, flow ceases and starts discontinuously with the 
closing or opening of the constrictions once adhesion is 
considered. 

Our present analysis does not include plasticity, because 
it will mostly affect the points in earliest contact, namely 
the peak heights where the local pressures are highest. 
We are mostly concerned with the areas with the most 
tenuous contact, and lowest contact pressures. Lateral 
plastic material flow from the peaks might slightly alter 
the detailed shape of the saddle point but such effects are 
beyond the scope of the present work. 

The most problematic approximations for the contin¬ 
uation of the power law down to the nanoscopic vicin¬ 
ity of the final constriction are most likely that we ne¬ 
glect to consider (i) the existence of a slip length (and the 
associated deviation from a Poisseuille flow profile), (ii) 
confinement-induced viscosity change (via the orientation 
of molecules), as well as (iii) the non-smoothness of sur¬ 
faces at nano-scales (where the continuum approximation 
breaks down). 

However, we note that shear-thinning does not come 
into play for the flow analysis, not even very close to the 
percolation threshold. This is because the maximum shear 
rate of a Poisseuille flow through a constriction is propor¬ 
tional to the gap height divided by the channel length. As 
a consequence, the maximum shear rate disappears with 
/C-x w ith ( — X ~ 0-6, which we see confirmed by our 
simulation results. 


Conclusions. — In this study we find that adjacent 
contact patches merge or break up discontinuously when 
finite-range interactions act in addition to hard-wall re¬ 
pulsion. This result has two main implications: First, 
in most cases, contact patches cannot be separated by 
a marginal distance or joined via arbitrarily thin ridges, 
which might explain why classical experiments [3] visual¬ 
izing the contact between two nominally flat solids as well 
as recent computer simulations including adhesion (25 


showed clearly separated contact patches, while large-scale 
simulations of non-adhesive solids [5jj9] usually necessitate 
excessively high resolution to determine the connectedness 
of contact patches. Second, the contact-patch coalesence 
and break-up occur through first-order instabilities imply¬ 
ing multistability. Thus, coalesence and break-up dynam¬ 
ics constitute a dissipation mechanism for Coulomb fric¬ 
tion as any other instability induced by the relative slid¬ 
ing motion of two solids 26 27 . Moreover, multistability 


entails thermal aging, which could be perceived as bulk 
relaxation and lead, for example, to a non-negligible time 
dependence of the leak rate of seals near percolation. 

In reference to Persson’s contact mechanics [2] and his 
single-junction leak-rate theory 14 , we also find two main 
implications. First, we demonstrate that the contact me¬ 
chanics of surfaces with single-wavelength roughness and 
their leakage rate are highly sensitive to the way in which 
local maxima and minima are arranged. For example, 
our triangular and hexagonal models seal at very different 
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normal loads and relative contact areas, although the two 
geometries have identical height spectra. It might be pos¬ 
sible to incorporate such effects — effectively representing 
a correlation of the phases for different h( q) at fixed |q| - 
into the theory via appropriate generalizations. Second, 
in the absence of adhesion, we find that the length and the 
width of a constriction are not similar in magnitude as as¬ 
sumed in Ref. 14 . Instead they disappear with different 


power laws as L c is approached. 

When adhesion is absent, the contact mechanics of 
“well-behaved” saddle points shows universal behavior 
near the percolation point. There, the gap is entirely de¬ 
termined by the sign of the reduced load and scaling fac¬ 
tors. Another aspect of the universal behavior is that the 
pressure profile on a newly-formed contact ridge along the 
transverse direction seems identical to that of a Hertzian 
contact. We therefore see the possibility for closed-form 
analytical solutions of the contact mechanics of a saddle 
point, even if formulating the boundary conditions could 
be difficult. 

Finally, we find for idealized conditions frequently as¬ 
sumed in continuum approaches (zero slip length, no con¬ 
finement effects on viscosity, hard-wall repulsion only) 
that the Reynolds flow through a critical constriction fol¬ 
lows (L C — L)P with [3 = 3.45. Since, for the given assump¬ 
tions, most fluid pressure drops near a single constriction 


in the immediate vicinity of the percolation point 14 , the 


leak rate of a seal that is pressed against a randomly rough 
macroscopic solid would depend on load with the same 
power law as the isolated constriction. Despite this knowl¬ 
edge it is quite difficult to predict the absolute current 
because the detailed topography of the last constriction 
has a tremendous effect on prefactors but is usually un¬ 
known. Still, the percolation exponent would result from 
local contact mechanics and not, as commonly assumed in 


the theory of percolation 28 , from the disorder at large 


length scales. Whether the classical view is reestablished 
once adhesion is included remains to be seen. 


We thank the Jiilich Supercomputing Centre for com¬ 
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